/*
 * \file CplexColorRD.cpp
 * \brief Implementación de la clase CplexColorRD
 *
 *  Created on: Jul 29, 2014
 *      Author: correa
 */

#include <CplexColorRD.h>
#include <Separator.h>
#include <SeparatorParameters.h>
#include <RDCliquer.h>
#include <HeuristicasGolosas.h>
#include <Reserva.h>
#include <TimeLimit.h>
#include <string>
#include <list>

//macro que se necesita para portabilidad, parece.
ILOSTLBEGIN

//macro para implementar fácilmente el callback para ramificación.
ILOBRANCHCALLBACK4(StabRepBranching, CplexRepBranching *, _branching, CplexRepSearchTree *, _tree, IloNumVarArray, _x, ArrayOfRepresentaciones&, _array)
{
	if (getNbranches() < 2)
		return;

	int bestr = -1;
	double bestvalue = 1.0;
	int bestv = -1;
	int r;
	for (int r = 0; r < _array.size(); r++) {
//	for (r = _array.size()-1; r >= 0; r--) {
		double max = 1.0-getValue(_x[_array[r][r]]);
		if (max > 0.55) {
			max = max/2.0;
			int v = 0;
			double xv = max;
			for (int i = r+1; i < _array.size(); ++i)
				if (_array[r][i] >= 0) {
					double nxv = fabs(getValue(_x[_array[r][i]])-max);
					if (nxv < xv) {
						v = i;
						xv = nxv;
					}
				}

			xv = xv/max;
			if (xv < bestvalue) {
				bestvalue = xv;
				bestr = r;
				bestv = v;
			}

			if (bestvalue < 0.1)
				break;
		}
		else if (max > 0.45 && max < 0.55)
		{
			bestv = -1;
			bestr = r;
			bestvalue = fabs(getValue(_x[_array[r][r]])-0.5);
			break;
		}
	}

	if (bestv < 0) {
		for (r--; r >= 0; r--) {
			double val = fabs(getValue(_x[_array[r][r]])-0.5);
			if (val < bestvalue) {
				bestr = r;
				bestvalue = val;
			}
		}

		cout << "bestr=" << bestr << " bestval=" << bestvalue << endl;
		makeBranch(_x[_array[bestr][bestr]], 0.0, IloCplex::BranchDown, getObjValue(), _tree->bajar(bestr, getNodeData()));
		makeBranch(_x[_array[bestr][bestr]], 1.0, IloCplex::BranchUp, getObjValue(), _tree->subir(bestr, getNodeData()));
		return;
	}

	if (1.0-getValue(_x[_array[bestr][bestr]]) - getValue(_x[_array[bestr][bestv]]) < 0.3) {
		cout << "bestr=" << bestr << " bestv=" << bestv << " bestval=" << bestvalue << endl;
		cout << " valuer=" << 1.0-getValue(_x[_array[bestr][bestr]]) << " valuev=" << getValue(_x[_array[bestr][bestv]]) << endl;
	}
	IloRangeArray constraints(getEnv());
	list<DesigualdadArray<RepRangeDesigualdad> *>& desigualdadesEqual = _branching->branchEqual(bestr, bestv);
	while (!desigualdadesEqual.empty())
	{
		CplexRepDesigualdad * d = (CplexRepDesigualdad *) desigualdadesEqual.back();
		desigualdadesEqual.pop_back();
		constraints.add(*d);
//		d->print();
		delete d;
	}
	IloCplex::MIPCallbackI::NodeData * data = _tree->identificar(bestr, bestv, getNodeData());
	_tree->dominados(bestr, data, constraints);
	_tree->dominados(bestv, data, constraints);
	makeBranch(constraints, getObjValue(), data);

	data = _tree->ligar(bestr, bestv, getNodeData());
	constraints.clear();

	bool dr = _tree->dominados(bestr, data, constraints);
	bool dv = _tree->dominados(bestv, data, constraints);
	if (!dr && !dv)
	{
		IloNumVarArray vars(getEnv());
		vars.add(_x[_array[bestr][bestr]]);
		vars.add(_x[_array[bestr][bestv]]);
		IloNumArray bounds(getEnv());
		bounds.add(2, 0.0);
		IloCplex::BranchDirectionArray dirs(getEnv());
		dirs.add(2, IloCplex::BranchDown);
		makeBranch(vars, bounds, dirs, getObjValue(), data);
	}
	else
	{
		list<DesigualdadArray<RepRangeDesigualdad> *>& desigualdadesDiffer = _branching->branchDiffer(bestr, bestv);
		while (!desigualdadesDiffer.empty())
		{
			CplexRepDesigualdad * d = (CplexRepDesigualdad *) desigualdadesDiffer.back();
			desigualdadesDiffer.pop_back();
			constraints.add(*d);
//			d->print();
			delete d;
		}
		makeBranch(constraints, getObjValue());
	}
}

//macro para implementar fácilmente el callback para generar cortes.
ILOUSERCUTCALLBACK5(RepCutGenerator, vector<CplexRepSeparator *>, _separator, CplexRepSearchTree *, _tree, IloNumVarArray, _x, ArrayOfRepresentaciones&, _array, RepCutGeneratorHistory, _history)
{
	CplexRepSolucion solucion(this, _x, _array);
	int n = _array.size();
	int nnodes = getNnodes();
	IloCplex::MIPCallbackI::NodeData * data = getNodeData();
	if (data == NULL)
	{
		setNodeData(_tree->root());
		data = _tree->root();
	}
	if (nnodes > 0)
		_tree->configGrafo(data);

	int s = 0;
	for (int r = 0; r < n && s < 2; r++) {
		Solucion& subsolucion = *solucion.newInstance(r, n);
		if (data != _history.lastNode && nnodes == 0)
			_separator[r]->restart();
		if (getValue(_x[_array[r][r]]) < 0.49 && subsolucion.nfracVertices() >= _array[r].size()/5) {
			if (getValue(_x[_array[r][r]]) > 0.9)
			{
				subsolucion.print();
				subsolucion.fracVertices().print();
				cout << "há " << subsolucion.nZeroes() << " zeros" << endl;
			}

			list<Desigualdad *> desigualdades;
			if (nnodes == 0 || nnodes % 5 == 0)
				desigualdades = _separator[r]->separar(subsolucion);
			else {
				//	getEnv().out() << "Solucion de valor " << getObjValue() << " y best bound " << getBestObjValue() << "..." << endl;
				//	getEnv().out() << "La solución es [";
				//	for (int i = 0; i < _x.getSize(); i++)
				//		getEnv().out() << " (" << i << ")" << solucion[i];
				//	getEnv().out() << " ]" << endl;

				desigualdades = _separator[r]->separarReserva(subsolucion);
			}
			if (!desigualdades.empty())
				s++;

			int nclique = 0;
			int nrango = 0;
			int ncorte = 0;
			while (!desigualdades.empty())
			{
				CplexSubRepDesigualdad * d = (CplexSubRepDesigualdad *) desigualdades.back();
				ncorte++;
				desigualdades.pop_back();
				if (d->isRank())
				{
					if (d->getRHS() == 1)
					{
						nclique++;
						addLocal(*d);
					}
					else
					{
						nrango++;
						if (nnodes == 0)
							add(*d, IloCplex::UseCutFilter);
						else
							addLocal(*d);
					}
				}
				else
					if (nnodes == 0)
						add(*d, IloCplex::UseCutFilter);
					else
						addLocal(*d);
				getEnv().out() << "-----> Corte con " << d->getSoporte().size() << " vertices: ";
				d->print();
				delete d;
			}
		}
		delete &subsolucion;
	}

	_history.lastNode = data;
}

Subproblema::Subproblema(IloEnv env, Grafo& g, VertexSet& hb, DiGrafo& o, ArrayOfRepresentaciones& array, int sz) : _grafo(g), _order(o), _array(array), _nh(0), _H(g.getN()), _Hbar(hb), _tamanho(hb.size()), _nreps(0), _rep(g.getN()), _auxset(g.getN()), _HS(g.getN()), _SH(g.getN()), _degree(new int[g.getN()]), _solucionEntera(env, sz)
{
	for (int i = 0; i < sz; ++i)
		_solucionEntera[i] = 0;
};

Subproblema::~Subproblema()
{
	 delete[] _degree;
}

bool Subproblema::isActiva(int v, IloNumVarArray& x)
{
	return (x[_array[v][v]].getLB() < 1.0e-05 && x[_array[v][v]].getUB() > 1.0-1.0e-05);
}

void Subproblema::activar(int v, IloNumVarArray& x) {
	_auxset.addAll();
	_auxset.removeAll(_grafo.vecinos(v));
	_auxset.remove(v);
	for (VertexSetIterator w = _auxset.begin(); *w < v && w != _auxset.end(); ++w)
		x[_array[*w][v]].setBounds(0, 1);
	x[_array[v][v]].setBounds(0, 1);
}

void Subproblema::desactivar(int v, IloNumVarArray& x) {
	_auxset.addAll();
	_auxset.removeAll(_grafo.vecinos(v));
	_auxset.remove(v);
	for (VertexSetIterator w = _auxset.begin(); *w < v && w != _auxset.end(); ++w)
		x[_array[*w][v]].setBounds(0, 0);
	x[_array[v][v]].setBounds(1, 1);
}

void Subproblema::activarTodos(IloNumVarArray& x, IloNumVarArray& y) {
	for (int w = 0; w < _H.capacity(); ++w)
	if (_H.contains(w))
	{
		activar(w, x);
		y[w].setBounds(0, 0);
	}
	else
	{
		desactivar(w, x);
		y[w].setBounds(0, 0);
	}
}

int Subproblema::buscarYDesactivar(IloNumVarArray& x, IloNumVarArray& y) {
	int ret = 0;
	VertexSet auxrep(_rep.capacity());
	VertexSet nbar(_H.capacity());

	for (VertexSetIterator w = _H.begin(); w != _H.end(); ++w)
	{
//		cout << "w=" << *w << " deg=" << _degree[*w] << " nreps=" << _nreps << " is act? " << x[_array[*w][*w]].getUB() << endl;

		if (_degree[*w] < _nreps)
		{
			ret++;
			if (isActiva(*w, x))
			{
				desactivar(*w, x);
				y[*w].setBounds(1, 1);

//				cout << "xww=" << operator [](_array[*w][*w]) << endl;

				if (_solucionEntera[_array[*w][*w]] <= 1.0e-05)
				{
					_solucionEntera[_array[*w][*w]] = 1;
					_rep.remove(*w);
					nbar.copyAll(_H);
					nbar.removeAll(_grafo.vecinos(*w));
					int rr = -1;
					for (VertexSetIterator v = nbar.begin(*w+1); v != nbar.end(); ++v)
					{
						if (_solucionEntera[_array[*w][*v]] >= 1.0-1.0e-05)
						{
							if (rr == -1)
							{
								rr = *v;
								_solucionEntera[_array[*v][*v]] = 0;
								_rep.add(rr);
							}
							else
								_solucionEntera[_array[rr][*v]] = 1;
							_solucionEntera[_array[*w][*v]] = 0;
						}
					}
				}
				else
				{
					auxrep.copyAll(_rep);
					auxrep.removeAll(_grafo.vecinos(*w));
					for (VertexSetIterator u = auxrep.begin(); *u < *w && u != auxrep.end(); ++u)
						_solucionEntera[_array[*u][*w]] = 0;
				}
			}
		}
	}

	return ret;
}

int Subproblema::assignx(IloCplex& cplex, IloNumVarArray& x, IloNumVarArray& y)
{
	cplex.getValues(x, _solucionEntera);

	_rep.removeAll();
	int p = _nreps;
	_nreps = 0;
	for (VertexSetIterator w = _H.begin(); w != _H.end(); ++w)
	{
		if (isActiva(*w,x) && _solucionEntera[_array[*w][*w]] <= 1.0e-5)
		{
			_rep.add(*w);
			_nreps++;
		}
	}

	if (_nreps >= p)
		buscarYDesactivar(x, y);

	return 0;
}

int Subproblema::addCandidato(IloRangeArray& constraints, IloNumVarArray& x, IloNumVarArray& y)
{
	int v = _Hbar.least();
	if (v >= 0)
		addCandidato(v, constraints, x, y);
	return v;
}

int Subproblema::addCandidato(int r, IloRangeArray& constraints, IloNumVarArray& x, IloNumVarArray& y)
{
	VertexSet auxrep(_rep.capacity());
	VertexSet nbar(_H.capacity());

	int act = 0;
	VertexSet auxH(_H);
	auxH.retainAll(_grafo.vecinos(r));
	_degree[r] = 0;
	for (VertexSetIterator w = auxH.begin(); w != auxH.end(); ++w)
	{
//		cout << "r=" << r << " w=" << *w << " degw=" << _degree[*w] << " reps=" << _nreps << endl;

		_degree[r]++;
		if (++_degree[*w] == _nreps && !isActiva(*w, x))
		{
			activar(*w, x);
			y[*w].setBounds(0, 0);
			act++;

			auxrep.copyAll(_rep);
			auxrep.removeAll(_grafo.vecinos(*w));
			VertexSetIterator u = auxrep.begin();
			if (u == auxrep.end()) // ocurre si sólo hay 1 rep porque todos estan desactivados
			{
				_solucionEntera[_array[*w][*w]] = 0;
				_rep.add(*w);
			}
			else for (; u != auxrep.end(); ++u)
			{
				nbar.copyAll(_H);
				nbar.removeAll(_grafo.vecinos(*u));
				bool urep = true;
				for (VertexSetIterator v = nbar.begin(*u+1); urep && v != nbar.end(); ++v)
					if (_solucionEntera[_array[*u][*v]] >= 1.0-1.0e-05 && _grafo.isArista(*w, *v))
					{
						urep = false;
						break;
					}
				if (urep)
					if (*u < *w)
					{
						_solucionEntera[_array[*u][*w]] = 1;
						_solucionEntera[_array[*w][*w]] = 1;
					}
					else
					{
						_solucionEntera[_array[*u][*u]] = 1;
						_rep.remove(*u);
						_solucionEntera[_array[*w][*w]] = 0;
						_rep.add(*w);
						_solucionEntera[_array[*w][*u]] = 1;
						for (VertexSetIterator v = nbar.begin(*u+1); v != nbar.end(); ++v)
							if (_solucionEntera[_array[*u][*v]] >= 1.0-1.0e-5)
							{
								_solucionEntera[_array[*u][*v]] = 0;
								_solucionEntera[_array[*w][*v]] = 1;
							}
					}
			}
		}
	}

	if (act > 0)
		greedify(constraints, x, y);

	_nh++;
	_H.add(r);
	_Hbar.remove(r);

	if (_degree[r] >= _nreps)
	{
		activar(r,x);
		y[r].setBounds(0, 0);
		_solucionEntera[_array[r][r]] = 1;
//		_rep.add(r);
//		_nreps++;
	}

	return 1;
}

int Subproblema::extendWithStab(IloRangeArray& constraints, IloNumVarArray& x, IloNumVarArray& y)
{
	int v = _Hbar.least();
	if (v < 0)
		return 0;
	return extendWithStab(v, constraints, x, y);
}

int Subproblema::extendWithStab(int r, IloRangeArray& constraints, IloNumVarArray& x, IloNumVarArray& y)
{
//	cout << "r=" << r << endl;

	int set[_grafo.getN()];
	int ns = 0;
	int rr = r;
	bool ract = isActiva(rr, x);
	set[ns++] = r;

	VertexSet auxH(_Hbar);//(_grafo.getN());
	auxH.removeAll(_grafo.vecinos(r));
	auxH.remove(r);

	int v;
	while ((v = auxH.least()) >= 0) {
//		cout << "v=" << v << endl;
		if (v < r)
			r = v;
		set[ns++] = v;

		auxH.removeAll(_grafo.vecinos(v));
		auxH.remove(v);
	}

	VertexSet auxrep(_rep.capacity());
	VertexSet nbar(_H.capacity());

	int act = 0;
	for (int i = 0; i < ns; ++i)
		if (set[i] != rr || !ract)
		{
			auxH.copyAll(_H);
			auxH.retainAll(_grafo.vecinos(set[i]));
			_degree[set[i]] = 0;
			for (VertexSetIterator w = auxH.begin(); w != auxH.end(); ++w)
			{
//				cout << "seti=" << set[i] << " w=" << *w << " degw=" << _degree[*w] << " reps=" << _nreps << endl;

				_degree[set[i]]++;
				if (++_degree[*w] == _nreps && !isActiva(*w, x))
				{
					activar(*w, x);
					y[*w].setBounds(0, 0);
					act++;

					auxrep.copyAll(_rep);
					auxrep.removeAll(_grafo.vecinos(*w));
					VertexSetIterator u = auxrep.begin();
					bool urep = false;
					if (u == auxrep.end()) // ocurre si sólo hay 1 rep porque todos estan desactivados
					{
						_solucionEntera[_array[*w][*w]] = 0;
						_rep.add(*w);
					}
					else for (; !urep && u != auxrep.end(); ++u)
					{
						nbar.copyAll(_H);
						nbar.removeAll(_grafo.vecinos(*u));
						urep = true;
						for (VertexSetIterator v = nbar.begin(*u+1); urep && v != nbar.end(); ++v)
							if (_solucionEntera[_array[*u][*v]] >= 1.0-1.0e-05 && _grafo.isArista(*w, *v))
								urep = false;
						if (urep)
							if (*u < *w)
							{
								_solucionEntera[_array[*u][*w]] = 1;
								_solucionEntera[_array[*w][*w]] = 1;
							}
							else
							{
								_solucionEntera[_array[*u][*u]] = 1;
								_rep.remove(*u);
								_solucionEntera[_array[*w][*w]] = 0;
								_rep.add(*w);
								_solucionEntera[_array[*w][*u]] = 1;
								for (VertexSetIterator v = nbar.begin(*u+1); v != nbar.end(); ++v)
									if (_solucionEntera[_array[*u][*v]] >= 1.0-1.0e-5)
									{
										_solucionEntera[_array[*u][*v]] = 0;
										_solucionEntera[_array[*w][*v]] = 1;
									}
							}
					}
				}
			}
		}

	for (int i = 0; i < ns; ++i)
	{
		_H.add(set[i]);
		_Hbar.remove(set[i]);

		if (_degree[set[i]] >= _nreps)
		{
			if (set[i] != rr || !ract)
			{
				activar(set[i],x);
				y[set[i]].setBounds(0, 0);
			}
			if (set[i] != r) {
//				cout << " ";
				_solucionEntera[_array[r][set[i]]] = 1;
				_solucionEntera[_array[set[i]][set[i]]] = 1;
			}
			else
			{
//				cout << " r";
				_solucionEntera[_array[set[i]][set[i]]] = 0;
			}
		}
//		else
//			cout << " ex";
//		cout << set[i] << "(deg=" << _degree[set[i]] << ")";
	}

	if (act > 0)
		greedify(constraints, x, y);

	if (ns > 0)
	{
		ns--;
		_nh+=ns;
		_rep.add(r);
//		 if (!ract)
			 _nreps++;
//		 else
//		 {
//			 if (r != rr)
//				 _rep.remove(rr);
//		 }
	}

//	cout << "\nacrescentados " << ns << " vertices com " << _nreps << " representantes." << endl;

	return ns;
}

void Subproblema::addChromoBoundConst(int objValue, IloRangeArray& constraints, IloNumVarArray& x, IloNumVarArray& y)
{
	IloRange cromconst(constraints.getEnv(), 0, objValue);
	for (VertexSetIterator itr = _H.begin(); itr != _H.end(); ++itr) {
		_auxset.copyCompl(_grafo.vecinos(*itr));
		_auxset.retainAll(_H);
		for (VertexSetIterator u = _auxset.begin(*itr+1); u != _auxset.end(); ++u)
		{
			cromconst.setLinearCoef(x[_array[*itr][*u]], 1);
			cout << " + " << x[_array[*itr][*u]].getName();
		}
		cromconst.setLinearCoef(y[*itr], 1);
		cout << " + " << y[*itr].getName();
	}
	cout << " <= " << objValue << endl;
	constraints.add(cromconst);
}

void Subproblema::addChromoConst(int objValue, IloRangeArray& constraints, IloNumVarArray& x)
{
	_HS.copyAll(_H);
	_SH.removeAll();
	int ns = 0;
	for (VertexSetIterator v = _H.begin(); v != _H.end(); ++v)
		if (_SH.isSubsetOf(_grafo.vecinos(*v)))
		{
			_SH.add(*v);
			ns++;
			_HS.remove(*v);
		}
	if (ns < _nreps) {
		IloRange colorconst(constraints.getEnv(), -objValue, _grafo.getN());

		for (VertexSetIterator v = _HS.begin(); v != _HS.end(); ++v) {
			colorconst.setLinearCoef(x[_array[*v][*v]], -1);
			_auxset.copyCompl(_H);
			_auxset.removeAll(_grafo.vecinos(*v));
			for (VertexSetIterator u = _auxset.begin(); *u < *v && u != _auxset.end(); ++u)
//				if (!_grafo.isArista(*u, *v))
					colorconst.setLinearCoef(x[_array[*u][*v]], 1);
		}
		constraints.add(colorconst);
	}

	// desigualdad de clique para SH
	if (ns > 2)
		addCliqueConst(_SH, constraints, x);
}

void Subproblema::addCliqueConst(VertexSet& clique, IloRangeArray& constraints, IloNumVarArray& x)
{
	for (VertexSetIterator r = _Hbar.begin(); r != _Hbar.end(); ++r) {
		_auxset.copyCompl(_grafo.vecinos(*r));
		_auxset.retainAll(clique);
		IloRange colorconst(constraints.getEnv(), 0, 1);
		colorconst.setLinearCoef(x[_array[*r][*r]], 1);
		for (VertexSetIterator v = _auxset.begin(*r+1); v != _auxset.end(); ++v)
			colorconst.setLinearCoef(x[_array[*r][*v]], 1);
		constraints.add(colorconst);
	}
}

void Subproblema::addNosymConst(int r1, int r2, int v, IloRangeArray& constraints, IloNumVarArray& x)
{
	IloRange nosymconst(constraints.getEnv(), 0, _grafo.getN());
	nosymconst.setLinearCoef(x[_array[r1][v]], -1);
	cout << "- " << x[_array[r1][v]].getName();
	nosymconst.setLinearCoef(x[_array[r2][r2]], 1);
	cout << " + " << x[_array[r2][r2]].getName();

	_auxset.copyCompl(_grafo.vecinos(r2));
	_auxset.retainAll(_grafo.vecinos(v));
	for (VertexSetIterator w = _auxset.begin(r2); w != _auxset.end(); ++w) {
		nosymconst.setLinearCoef(x[_array[r2][*w]], 1);
		cout << " + " << x[_array[r2][*w]].getName();
	}
	cout << " >= " << nosymconst.getLB() << endl;
	constraints.add(nosymconst);
}

void Subproblema::addNosymConsts(int r1, IloRangeArray& constraints, IloNumVarArray& x)
{
	VertexSet novecinos(_grafo.getN());
	for (VertexSetIterator r2 = _rep.begin(); r2 != _rep.end(); ++r2)
		if (*r2 != r1 && !_order.isArco(*r2,r1))
		{
			_order.addArco(r1,*r2);
			novecinos.copyCompl(_grafo.vecinos(r1));
			novecinos.removeAll(_grafo.vecinos(*r2));
			novecinos.retainAll(_H);
			novecinos.remove(r1);
			novecinos.remove(*r2);
			int max = r1 > *r2 ? r1 : *r2;
			for (VertexSetIterator v = novecinos.begin(max); v != novecinos.end(); ++v)
				addNosymConst(r1, *r2, *v, constraints, x);
		}
}

void Subproblema::createNosingletonConst(int v, IloRange& constraint, IloNumVarArray& x)
{
	_auxset.copyAll(_H);
	_auxset.removeAll(_grafo.vecinos(v));
	constraint.setLB(1);
	constraint.setUB(_grafo.getN());
	constraint.setLinearCoef(x[_array[v][v]], 1);
	for (VertexSetIterator w = _auxset.begin(v+1); w != _auxset.end(); ++w)
		constraint.setLinearCoef(x[_array[v][*w]], 1);
}

IloNumArray& Subproblema::solucionEntera()
{
	return _solucionEntera;
}

void Subproblema::greedify(IloRangeArray& constraints, IloNumVarArray& x, IloNumVarArray& y) {
	Lista& representantes = _rep.lista();
//	cout << "representantes: ";
//	representantes.print();

	VertexSet novecinos1(_grafo.getN());
	VertexSet novecinos2(_grafo.getN());
	VertexSet moved(_grafo.getN());

	// greedify vertices already in _H
	for (int ir1 = 0; ir1 < representantes.cantidad(); ir1++) {
		int r1 = representantes[ir1];
		novecinos1.copyCompl(_grafo.vecinos(r1));
		novecinos1.retainAll(_H);
		novecinos1.removeAll(moved);
		novecinos1.remove(r1);
		for (VertexSetIterator v = novecinos1.begin(r1); v != novecinos1.end(); ++v)
			if (_solucionEntera[_array[r1][*v]] >1.0-1.0e-05) {
				bool hayvecino = true;
				int r2 = r1;
				for (int ir2 = 0; hayvecino && ir2 < representantes.cantidad(); ir2++) {
					int r = representantes[ir2];
					if (r < *v && _order.vecinos(r2).contains(r) && _order.vecinos(r1).contains(r) && !_grafo.isArista(r, *v))
					{
						novecinos2.copyCompl(_grafo.vecinos(r2));
						novecinos2.retainAll(_grafo.vecinos(*v));
						novecinos2.remove(r2);
						novecinos2.retainAll(_H);
						VertexSetIterator w = novecinos2.begin(r2);
						hayvecino = (w == novecinos2.end());
						for (; !hayvecino && w != novecinos2.end(); ++w)
							hayvecino = (_solucionEntera[_array[r2][*w]] >1.0-1.0e-05);
					}
					if (!hayvecino)
						r2 = r;
				}
				if (!hayvecino) {
					moved.add(*v);
					_solucionEntera[_array[r1][*v]] = 0;
					_solucionEntera[_array[r2][*v]] = 1;

					addNosymConst(r1, r2, *v, constraints, x);
				}
			}
	}
}

void Subproblema::print(IloNumVarArray& y)
{
	cout << "================================" << endl;
	cout << "H (tam=" << _H.size() << "): ";
	_H.print();
	cout << "Reps (tam=" << _rep.size() << "): ";
	_rep.print();
//	for (VertexSetIterator u = _rep.begin(); u != _rep.end(); ++u)
//	{
//		cout << *u << ":";
//		_auxset.copyAll(_H);
//		_auxset.removeAll(_grafo.vecinos(*u));
//		for (VertexSetIterator v = _auxset.begin(*u); v != _auxset.end(); ++v)
//			if (_solucionEntera[_array[*u][*v]] >= 1.0-1.0e-05)
//				cout << " " << *v;
//		cout << endl;
//	}
	cout << "Inactivos: [";
	for (VertexSetIterator u = _H.begin(); u != _H.end(); ++u)
		if (y[*u].getUB() >= 1.0-1.0e-05)
			cout << " " << *u;
	cout << " ]" << endl;
	cout << "================================" << endl;
}

CplexStabColorModel::CplexStabColorModel(IloModel model, Grafo& _grafo, const int _k, MatrixOfRepresentaciones& _mapa, ArrayOfRepresentaciones& _array, SeparatorParameters& parameters) : IloCplex(model), _n(_grafo.getN()),_x(getEnv()), _y(getEnv()), _separator(_grafo.getN())
{
	//variable definition
	for (int r = 0; r < _n; ++r)
		for (int u = 0; u < _mapa[r].size(); u++)
			_x.add(IloIntVar(getEnv(), 0, 1));
	// variables describing vertex inactivity
	for (int r = 0; r < _n; ++r)
		_y.add(IloIntVar(getEnv(), 0, 1));

	//objective function
	IloNumExpr fobj(getEnv());
	for (int r = 0; r < _n; ++r) {
		string name("x"+to_string(r)+"_"+to_string(r));
		_x[_array[r][r]].setName(name.c_str());
		_x[_array[r][r]].setBounds(1, 1);
		for (int u = r+1; u < _n; ++u)
			if (!_grafo.isArista(r, u)) {
				string name("_x"+to_string(r)+"_"+to_string(u));
				_x[_array[r][u]].setName(name.c_str());
				_x[_array[r][u]].setBounds(0, 1);
				fobj += _x[_array[r][u]];
			}
		string yname("y"+to_string(r));
		_y[r].setName(yname.c_str());
		_y[r].setBounds(0, 0);
		fobj += _y[r];
	}
	getModel().add(IloMaximize(getEnv(), fobj));

	//seteamos los parámetros
	//		cplex.setParam(IloCplex::NodeLim, 0);
	setParam(IloCplex::Cliques, -1);
	setParam(IloCplex::ZeroHalfCuts, -1);
	setParam(IloCplex::FracCuts, -1);
	setParam(IloCplex::ImplBd, -1);
	//		cplex.setParam(IloCplex::TiLim,TimeLimit::MAX_SECONDS);
	//cplex.setParam(IloCplex.DoubleParam.TiLim, 60.0);

	//		cplex.setParam(IloCplex::Probe, -1);
	//		cplex.setParam(IloCplex::MIPEmphasis, IloCplex::BestBound);
	//		cplex.setParam(IloCplex::BndStrenInd, 0);
	//		cplex.setParam(IloCplex::RepeatPresolve, 0);
	//		cplex.setParam(IloCplex::HeurFreq, -1);
	//		cplex.setParam(IloCplex::FPHeur, -1);
	//		cplex.setParam(IloCplex::NodeSel, IloCplex::DFS);

	setParam(IloCplex::MIPDisplay, 2);
	setParam(IloCplex::Threads, 1);

	//edge/clique constraints
	StabBranching * rbranching[_n];

	IloRangeArray constraints(getEnv());
	VertexSet mask(_n);
	VertexSet vset(_n);
	for (int r = 0; r < _n; ++r) {
		vset.copyCompl(_grafo.vecinos(r));
		mask.add(r);
		vset.removeAll(mask);

		rbranching[r] = new CplexSubRepBranching(_grafo, r, _x, _mapa, _array, vset);

		list<Desigualdad *> desigualdades = rbranching[r]->branchNone();
		while (!desigualdades.empty())
		{
			CplexSubRepDesigualdad * d = (CplexSubRepDesigualdad *) desigualdades.back();
			//				d->print();
			desigualdades.pop_back();
			constraints.add(*d);
			delete d;
		}

		//			VertexSet aristas(n);
		//			for (VertexSetIterator v = novecinos.begin(); v != novecinos.end(); ++v) {
		//				aristas.copyAll(novecinos);
		//				aristas.retainAll(_grafo.vecinos(*v));
		//				VertexSetIterator u = aristas.begin();
		//				if (u == aristas.end())
		//					constraints.add(x[_array[r][*v]] + x[_array[r][r]] <= 1);
		//				else for (; u != aristas.end() && *u < *v; ++u)
		//					constraints.add(x[_array[r][*u]] + x[_array[r][*v]] + x[_array[r][r]] <= 1);
		//			}

		// add cover constraints
		IloRange coverconst(getEnv(), -1, 0);
		coverconst.setLinearCoef(_x[_array[r][r]], -1);

		vset.copyCompl(_grafo.vecinos(r));
		for (VertexSetIterator v = vset.begin(); v != vset.end() && *v < r; ++v)
			coverconst.setLinearCoef(_x[_array[*v][r]], 1);
		constraints.add(coverconst);

		getModel().add(constraints);
	}

	if (_k < _n)
	{
		constraints.clear();
		IloRange kconst(getEnv(), _n-_k, _n);
		for (int r = 0; r < _n; ++r)
			kconst.setLinearCoef(_x[_array[r][r]], 1);
		constraints.add(kconst);
		getModel().add(constraints);
	}

	// callback de generación de cortes y de _branching
	for (int r = 0; r < _n; ++r)
		_separator[r] = new CplexRepSeparator(_grafo, r, _x, _mapa, _array, parameters);
	_tree = new CplexRepSearchTree(_grafo, _x, _mapa, _array);
	_branching = new CplexRepBranching(_grafo, *_tree);
	use(RepCutGenerator(getEnv(), _separator, _tree, _x, _array, { NULL }));
	use(StabRepBranching(getEnv(), _branching, _tree, _x, _array));
};

CplexStabColorModel::~CplexStabColorModel()
{
	for (int r = 0; r < _n; ++r)
		delete _separator[r];
	delete _tree;
	delete _branching;
}

DoubleValueArray* CplexColorRD::solve(SeparatorParameters& parameters)
{
	DoubleValueArray * ret = 0;

	try
	{
		IloEnv env;
		IloModel model(env);
		CplexStabColorModel cplex(model, _grafo, _k, _mapa, _array, parameters);

		int n = _grafo.getN();
		VertexSet vset(n);
		IloRangeArray constraints(env);
		IloNumVarArray& x = cplex.getX();
		IloNumVarArray& y = cplex.getY();

		// russian dolls method
		DiGrafo order(_grafo.getN());
		order.addAllVertices();

		// clique maxima con cobertura por conjuntos independientes
		vset.addAll();
		VertexSet clique(n);

		int cm;
		clique.removeAll();
		cm = cliqueMaximal(&vset, &clique);
		cout << "clique maximal: ";
		clique.print();

		if (cm < 2)
			return ret;

		Subproblema subproblem(env, _grafo, vset, order, _array, x.getSize());

		constraints.clear();

		for (VertexSetIterator itv = clique.begin(); subproblem.ncolores() < _k && itv != clique.end(); ++itv) {
			subproblem.addCandidato(*itv, constraints, x, y);
			subproblem.extendWithStab(*itv, constraints, x, y);
		}

		if (subproblem.tamanoMuneca() < subproblem.tamano())
		{
			// desigualdad de clique para la clique maxima
			if (cm > 2)
				subproblem.addCliqueConst(clique, constraints, x);
			model.add(constraints);

			vset.removeAll(clique);
		}
		else
			return ret;

		int mipStart;
		subproblem.activarTodos(x, y);

		subproblem.print(y);

		subproblem.buscarYDesactivar(x, y);
		int objValue = subproblem.tamanoMuneca()-subproblem.ncolores();

		subproblem.print(y);

		while (subproblem.tamanoMuneca() < subproblem.tamano()) {
			constraints.clear();
			int newv = subproblem.addCandidato(constraints, x, y);

			// desigualdad que impone numero cromatico
			subproblem.addChromoBoundConst(objValue+1, constraints, x, y);
			model.add(constraints);

			subproblem.print(y);

			// Optimize the problem and obtain solution.
			if (subproblem.ncolores() >= 2) {
				int prevObjValue = objValue;

				IloRange nosingconst(env, 0, 0);
				subproblem.createNosingletonConst(newv, nosingconst, x);
				model.add(nosingconst);

				constraints.clear();
				subproblem.addNosymConsts(newv, constraints, x);
				model.add(constraints);

//				for (int r = 0; r < n; ++r)
//					cplex.separator(r)->getReservaCliques()->haceActiva();

//				cplex.exportModel("color.lp");

				if (cplex.getNMIPStarts() == 0)
					mipStart = cplex.addMIPStart(x, subproblem.solucionEntera());
				else
					cplex.changeMIPStart(mipStart, x, subproblem.solucionEntera());
				int status = cplex.solve();
				IloAlgorithm::Status probstat = cplex.getStatus();
				if (!status && (probstat != IloAlgorithm::Status::Infeasible))
				{
					env.error() << "No se pudo resolver el modelo. " << status << endl;
					throw(-1);
				}
				env.out() << "Solution status = " << probstat << endl;

				if (probstat == IloAlgorithm::Status::Optimal)
					objValue = cplex.getObjValue();

				// greedify
				env.out() << "Objective value  = " << objValue << endl;
				if (prevObjValue == objValue && subproblem.ncolores() < _k)
				{
					constraints.clear();

					int agregados = subproblem.extendWithStab(newv, constraints, x, y);
					objValue += agregados;
					cout << "agregados con extensión (newv=" << newv << "): " << agregados << endl;
					subproblem.buscarYDesactivar(x, y);

					if (subproblem.tamanoMuneca() < subproblem.tamano())
					{
						// desigualdad que impone numero cromatico
						subproblem.addChromoBoundConst(objValue, constraints, x, y);
						subproblem.addChromoConst(objValue, constraints, x);
					}
					model.add(constraints);
				}
				else
				{
					int agregados = subproblem.assignx(cplex, x, y);
					objValue += agregados;
					cout << "agregados con recolor: " << agregados << endl;
				}
				env.out() << "Solution value  = " << subproblem.tamanoMuneca() << endl;
				env.out() << "Colors  = " << subproblem.ncolores() << endl;

				model.remove(nosingconst);
				subproblem.print(y);

				//					if (probstat == IloAlgorithm::Status::Infeasible)
				//						exit(0);
			}
		}

		env.end();
	} catch (IloAlgorithm::NotExtractedException& e)
	{
		cerr << "Not extractated exception caught: " << e << endl;
		cerr << e.getExtractable().getName() << endl;
	} catch (IloException& e)
	{
		cerr << "Concert exception caught: " << e << endl;
	} catch (...)
	{
		cerr << "Unknown exception caught" << endl;
	}

	return ret;
}

DoubleValueArray* CplexColorRD::solveLinearRelaxation()
{
	DoubleValueArray * ret = 0;
	IloEnv env;
	try
	{
		IloModel model(env);

		IloNumVarArray var(env);
		IloRangeArray con(env);

		int n = _grafo.getN();

		IloEnv env = model.getEnv();

		//variable definition
		IloNumVarArray & x = *(new IloNumVarArray(env));
		for (int r = 0; r < n; ++r) {
			x.add(IloNumVar(env, 0, 1));				// _x^r_r
			for (int u = 0; u < _mapa[r].size(); u++)
				x.add(IloNumVar(env, 0, 1));			// _x^r_u
		}

		//objective function
		IloNumExpr fobj(env);

		//constraints
		IloRangeArray constraints(env);
		VertexSet mask(n);
		VertexSet novecinos(n);
		for (int r = 0; r < n; ++r) {
			IloRange coverconst(env, 1, 1);
			fobj += x[_array[r][r]];
			coverconst.setLinearCoef(x[_array[r][r]], 1);

			novecinos.copyCompl(_grafo.vecinos(r));

			for (VertexSetIterator v = novecinos.begin(); v != novecinos.end() && *v < r; ++v)
				coverconst.setLinearCoef(x[_array[*v][r]], 1);

			mask.add(r);
			novecinos.removeAll(mask);
			VertexSet aristas(n);
			for (VertexSetIterator v = novecinos.begin(); v != novecinos.end(); ++v) {
				aristas.copyAll(novecinos);
				aristas.retainAll(_grafo.vecinos(*v));
				VertexSetIterator u = aristas.begin();
				if (u == aristas.end())
					constraints.add(x[_array[r][*v]] + x[_array[r][r]] <= 1);
				else for (; u != aristas.end() && *u < *v; ++u)
					constraints.add(x[_array[r][*u]] + x[_array[r][*v]] + x[_array[r][r]] <= 1);
			}
			constraints.add(coverconst);
		}
		model.add(constraints);
		model.add(IloMinimize(env, fobj));

		IloCplex cplex(model);

		// Optimize the problem and obtain solution.
		if (!cplex.solve())
		{
			env.error() << "No se pudo resolver el modelo." << endl;
			throw(-1);
		}

		env.out() << "Valor relajación lineal  = " << cplex.getObjValue() << endl;

		for (int i = 0; i < n; i++)
			env.out() << "_x[" << i << "]=" << cplex.getValue(x[i]) << endl;

	} catch (IloException& e)
	{
		cerr << "Concert exception caught: " << e << endl;
	} catch (...)
	{
		cerr << "Unknown exception caught" << endl;
	}

	env.end();

	return ret;
}
